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It is well known experimentally that well-quenched amorphous solids exhibit a plastic instability 
in the form of a catastrophic shear localization at a well defined value of the external strain. The 
instability may develop to a shear-band that in some cases is followed by a fracture. It is also known 
that the values of the yield-strain (and yield-stress), as well as the direction of the shear band 
with respect to the principal stress axis, vary considerably with variations in the external loading 
conditions. In this paper we present a microscopic theory of these phenomena for 2-dimensional 
athermal amorphous solids that are strained quasi-statically. We present analytic formulae for the 
yield-strains for different loading conditions, and well as for the angles of the shear bands. We 
explain that the external loading conditions determine the eigenvalues of the quadrupolar Eshelby 
inclusions which model the non-afnne displacement field. These inclusions model elementary plastic 
events and determine both the yield-strain and the direction of the shear-band. We show that the 
angles of the shear bands with respect to the principal stress axis are limited theoretically between 
30° and 60°. Available experimental data conform to this prediction. 
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I. INTRODUCTION 

The physics of amorphous solids raises some problems 
that do not exist in perfect crystals. While it is well 
known that disorder brings about (Anderson) localiza- 
tion of cigenfunctions |l|, it is commonly assumed that 
such localization is limited to eigenfunctions associated 
with high energies. Low-frequency eigenfunctions are be- 
lieved to be spatially extended. This picture fails in 
amorphous solids that are subjected to an external strain. 
Here one finds that the lowest energy cigenfunction can 
become localized. This localization can take the form 
of a plastic event that changes the local organization of 
particles through the creation of essential non-afnne dis- 
placement. Sometimes the plastic instability can exhibit 
itself as a system spanning event of shear localization 
along a line in 2-dimensional systems [2| and a plane in 
3-dimcnsional systems [3[ . These events are crucially im- 
portant in limiting the toughness of technologically im- 
portant materials like metallic glasses. These instabilities 
are the subject of this paper. 

To be as precise as possible, we will deal in this paper 
with 2-dimensional systems (although the 3-dimensional 
extension is available |3j) which are athermal (at temper- 
ature T = 0) and strained quasi-statically Q. We will 
consider systems which arc good glass formers, contain- 
ing N particles in a volume V, interacting via generic in- 
teraction potentials that are sufficiently smooth (with at 
least two continuous derivatives). The total energy can 
be written then in terms of the positions r\,r2, ■ ■ ■ r^ 
of these particles, U = U(ri,r2 ■ ■ ■ ,r?j). The Hessian 
matrix is defined as the second derivative |5f 



H: 



d 2 U(r 1 ,r 2 --- ,r N ) 
drjdrj 



(1) 



The Hessian is real and symmetric, and therefore can be 
diagonalized. Excluding Goldstonc modes whose eigen- 



values are zero due to continuous symmetries, all the 
other eigenvalues are real and positive as long as the sys- 
tem is mechanically stable. In equilibrium, without any 
mechanical loading, the eigenfunctions associated with 
the large eigenvalues are localized due to Anderson lo- 
calization, as it was mentioned above. However, all the 
cigenfunctions associated with the low eigenvalues (in- 
cluding all the excess modes that are typical to amor- 
phous solids) are spatially extended. It had been a major 
discovery of the last decade that at small values of the ex- 
ternal strain there appear "fundamental plastic events" 
in which the eigenvalues of Hessian matrix hit zero via 
a saddle node bifurcation, and simultaneously the as- 
sociated cigenfunctions get localized Q. This is a new 
mechanism for localization, different from the well known 
Anderson mechanism, and it is the basis for the under- 
standing of plastic instabilities in amorphous solids. It 
had been shown that at low values of the external stains 
the spatial ranges of these plastic events are system-size 
independent [7|, involving a small number of particles, 
of the order of 100. The displacement fields associated 
with fundamental plastic events are reliably modeled by 
Eshelby inclusions (see below for details). It had been a 
more recent discovery that at higher values of the exter- 
nal strains these fundamental plastic events can organize 
in highly correlated lines in 2-dimensional system or in 
planes in 3 dimensional systems [2, 3]. These highly cor- 
related densities of Eshelby inclusions are the microscopic 
manifestation of the shear localization. 

The key idea of the recent work on the shear local- 
ization is that at T — with quasi-static loading the 
preferred spatial organization of the density of Eshelby 
inclusions could be found by minimizing their total en- 
ergy in the strained system. In Refs. j^bl the theory was 
worked out for simple external shear in both 2 and 3 di- 
mensions. It was found that in the case of a pure external 
shear the shear-localization is realized by organizing the 



inclusions on a line (plane) in 2 (3) dimensions that are 
precisely at 45° with respect to the principal stress axis. 
An additional important result is an analytic prediction 
for the yield-strain (where shear localization occurs for 
the first time) in terms of the properties of the Eshelby 
solution for the fundamental plastic event (a single inclu- 
sion). 

It is well known however that for other loading condi- 
tions, e.g. uniaxial tension or compression, the value of 
the yield strain as well as the direction of the shear band 
can vary considerably [8[, indicating that the pure shear 
is a special case. It was reported recently in [9| that the 
difference can be related to the more general form of the 
Eshelby inclusion which models the fundamental plastic 
instability at different loading conditions. The aim of the 
present paper is to offer the theory in sufficient detail and 
to work out analytic predictions of the yield strain. The 
calculations involved are straightforward but sometimes 
cumbersome, and we make an effort to offer them in an 
optimally didactic way in the sequel. 

In Sect. [IT] we present results for AQS numerical sim- 
ulations in uniaxial compression and extension to pro- 
vide us with data on yield strains and directions of shear 
bands. We find the well known asymmetry in both mea- 
sures. In Sect. IHII we begin to develop the general theory 
for 2-dimcnsional Eshelby inclusions that is appropriate 
for the most general loading conditions. In Sect. IIVI we 
compute the all-important energy of N inclusions, and 
minimize it to find the preferred organizations and the 
resulting angle of the shear band. On the basis of this 
calculation we offer in Sect. |V]an analytic prediction for 
the yield strain. The final section IVII provides a short 
summary and a discussion. 

The reader should note that throughout this paper we 
will reserve the Greek letters a, /3, 7, . . . for tensor sub- 
scripts, while the Latin letters i,j, k, . . . will be used for 
particles and Eshelby inclusions. 



II. NUMERICAL SIMULATIONS 

To prepare data for the present analysis we have per- 
formed two-dimensional (2D) Molecular Dynamics simu- 
lations on a binary system which is known to be a good 
glas s former and has a quasi-crystalline ground state 
jlfj . Hl| . Each atom in the system is labeled as either 
"small" (S) or "large" (L) and all the particles interact via 
Lennard- Jones (LJ) potential. All distances \n — Vj\ are 
normalized by Xsl, the distance at which the LJ poten- 
tial between the two species becomes zero and the energy 
is normalized by csl which is the interaction energy be- 
tween the two species. Temperature was measured in 
units of esi/fcs where ks is Boltzmann's constant. For 
detailed information on the model potential and its prop- 
erties, we refer the reader to Ref [10|- The number of 
particles in our simulations is 10000 at a number density 
n = 0.985 with a particle ratio N L /N S = (I + V5)/4. 
The mode coupling temperature Tmct for this system is 



known to be close to 0.325. The mass of all particles is 

mo = I and time is normalized to to = \l esL^SL /mo- 
For the sake of computational efficiency, the interaction 
potential is smoothly truncated to zero along with its first 
two derivatives at a cut-off distance r c = 2.5. To prepare 
the glasses, we first start from a well equilibrated liquid 
at a high temperature of T = 1.2 which is supercooled to 
T = 0.35 at relatively fast quenching rate of 3.4 x 10~ 3 . 
Then, we equilibrated these supercooled liquids for times 
greater than 20T re i, where r ro i is the time taken for the 
self intermediate scattering function to approach 1% of 
its initial value. Lastly, following this equilibration, we 
quenched these supercooled liquids deep into the glassy 
regime at a temperature of T = 0.01 at a reduced quench 
rate of 3.2 x 10 -6 . Following this quench we took the 
glass to mechanical equilibrium (nearest energy minima) 
by a conjugate gradient energy minimization. After that 
we start loading our glasses under athermal quasi-static 
conditions. Two different loading protocols were applied, 
one is uniaxial compression and the other - uniaxial ex- 
tension. The simulations were followed until the yield 
strain j Y and slightly above to ascertain the angle of the 
shear band formed. The shear band is marked on the 
sample by coloring particles using their values of -D^ in 
[12| . We use a binary coloring scheme which means that 
when a given particle i has D 2 min [i) > A| L , it is colored 
black, else it is colored white. 

We note the asymmetry in j Y , close to 5.5% for com- 
pression and 3.5% for extension, and the asymmetry in 
angles, 46° and 54° respectively. The theory outlined in 
this paper provides analytic formulae for both measures, 
cf Eqs. CD) and JTSJ). 



III. THEORY OF THE FUNDAMENTAL 
PLASTIC EVENT 

In the past, phenomenological models were applied to 
this problem |13l - [l7| . but a microscopic approach was 
lacking. In this section we present the theory of 2- 
dimensional Eshelby inclusions for the general loading 
conditions. 



A. Two Dimensional Circular Inclusion 

Consider an elastic solid having a volume V and 
surface area S [Fig. [3J. The material will be assumed 
to be homogeneous with an elastic stiffness tensor given 
by Cap-ys- Let a sub- volume Vq with surface area Sq 
undergo a uniform permanent (inelastic) deformation, 
such as a structural phase transformation. The material 
inside Vq is referred to as an inclusion and the material 
outside is called the matrix. If we could remove this 
inclusion from its surrounding material then it would 
attain a state of a uniform strain and zero stress. Such 
a stress free strain is referred to as the eigen-strain e*o. 
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FIG. 2. Color Online: A typical stress strain curve obtained 
under AQS conditions for uniaxial tension (lower curve: red 
filled circles) and compression (upper curve: green crosses). 
The nonlinear response is clearly asymmetric with yield strain 
7 Y being sa 5.5% for compression and 3.5% for extension re- 
spectively. In Sect. [V] we provide an analytic formula for 
yield strain which predicts this asymmetry to a very high de- 
gree of accuracy. 



FIG. 1. Color Online: The shear band that occurs in a 2- 
dimensional amorphous solid upon uniaxial compression (a) 
and extension (b). The angle with respect to the principal 
stress component measured in compression is 46°±1°, whereas 
in the extension it is 54° ±1°. 



The eigen stress is then given by a* a o = Cap^seZg- 

In reality, the inclusion is surrounded by the matrix. 
Therefore, it is not able to reach the state of zero stress. 
Instead, both the inclusion and the matrix will deform 
and experience an elastic stress field. The Eshelby's 
transformed inclusion problem [18[ is to solve the stress, 
strain and displacement fields both in the inclusion and 
in the matrix. 

We consider a 2D circular inclusion that has been 
strained into an ellipse using an eigen-strain e* a and 
which allows for a volume change (e*„ 7^ 0). A general 
expression for such a tensor can be written in terms of 
a unit eigenvectors (n, fc) and corresponding eigenvalues 
(C„,Cfc) as follows: 



1/} = (nhahp + (kk a k/3 



(2) 




FIG. 3. Cartoon showing an elastic medium of volume V and 
surface area S. Inside the medium a small ellipsoidal region 
(volume Vo and surface area So) undergoes an irreversible 
(plastic) deformation. The material inside Vq is called as the 
inclusion and the material outside is referred to as the matrix. 



(Cn-Ck) 



t a/3 



(2h a h l 3 -8 a p) + 



(Cn + Cfc) 



? a f3 



*,0 

£ a/3 ' 



E a/3 



(3) 



Using orthogonality of the eigen directions: h a np 
k a kp = 6 a p, we get: 



where the traceless part eJg and the nonzero trace part 



* T 
6 a/3 are g lven as 



*.0 \^n Sfe I /_ - 5 x \ 

£ afJ = o (2«o"-/3 - <W) 



(Cn+Cfe) 



~a£ 



3q/3 



(4) 



We also assume that the system is acted upon by a 
homogeneous strain e^ that acts globally (which in 
our case also triggers the local transformation of the 
inclusion). This strained ellipsoidal inclusion feels a 
traction exerted by the surrounding clastic medium 
resulting in a constrained strain t c a a in the inclusion and 
also exerts a traction at the inclusion-matrix interface 
resulting in the originally unstrained surroundings 
developing a constrained strain field e c a a{X). 

The eigen-strain e*o in the inclusion is related to the 
constrained strain e c a p via the fourth order Eshelby tensor 

Safins'- 



-a0 






76 



(5) 



Now for an inclusion of arbitrary shape the constrained 
strain e^n, both the stress a c a o and the displacement field 
u c (X) inside the inclusion are in general functions of 
space. For ellipsoidal inclusions, however, it was shown 
by Eshelby [18|, [jjj that the Eshelby tensor and the con- 
strained stress and strain fields inside the inclusion be- 
come independent of space. We work here with a circular 
inclusion which is a special case of an ellipse and hence 
for such an inclusion, the Eshelby tensor S a p 1 s reads [20( 



4(A + 2 M ) 



4(A + 2//) 



(6) 

where A, [i are Lame coefficients. Plugging Eq ^ in Eq 
©, we get 



(A + 3fi) 
" aP ~ (A + 2 M ) C ^ T 2(A + 2/x)' 



(\ + V)_ *,T 






(7) 



The total stress, strain and displacement field inside 
the circular inclusion are then given by 



-a/3 



E a/3 T t a p 



I C _* I _oo /"-» / C ,,* I oo\ 

CT a/3 — CT a,3 ~~ °a,S + V = °a/37<H e 7<5 ~~ e -y<5 + e i$> 



C i U 



( e afj + e ^) X l3 



(8) 



where the superscript J indicates the inclusion. The eigen 
stress <7*g is related to the eigen strain e* « as 



7 a/3 



n * 

2 M e a/3 + ^™<W 



(9) 



where we have used the following definition of the 
fourth order elastic stiffness tensor C a p 1 s f° r an isotropic 
elastic medium: 



C a f3~fS = ^$af3S~/S + ^{SaySpS + SaSdp-,) 



(10) 



The stress in the inclusion can now be written down 
in terms of the independent variables using equations ((9|) 
as 



'a/3 



2/x 



-a0~ t a0~ rt a0 



+ A 



S a fi (11) 



B. Constrained Fields in the Matrix 

In the surrounding elastic matrix, the stress, strain and 
displacement fields are all explicit functions of space and 
can be written as 



c a/3 



C/3P0 = e a/3 (X) - 

u^{X) = ul fj {X)+u^ (12) 

The displacement field u c a {X) in the isotropic elastic 
medium will satisfy the Lame-Navier equation (without 
any body forces) [2l[ 



(P + A) 



d 2 !!*: 



d 2 u° 



/'- 



dX a dX 7 ^dXpdXp 







(13) 



The constrained fields in the inclusion will supply the 
boundary conditions for the displacement field in the 
matrix at the inclusion boundary. Also as r — > 00, the 
constrained displacement field will vanish. 



All solutions of Eq. 
bi-harmonic equation 



JU|) also obey the higher order 



dXsdXfsdXxdX) 



= V 4 i£ = 



VpOAp 



(14) 



Thus our objective is to construct from the radial so- 
lutions of the bi-Laplacian equation Eq. (|14p derivatives 
which also satisfy Eq. (p~3|) . Note that the bi-Laplacian 
equation is only a necessary (but not a sufficient) con- 
dition for the solutions and Eq. (JTHJ) still needs to be 
satisfied. 



C. Solution of the Lame-Navier Equation 

From the foregoing section, we note that the con- 
strained displacement field due to the Eshelby solution 



3 given as 




u a = e aP X P 


= 


' (A + /i) *, T (A + 3/x) +i0 " 
_(A + 2 M ) C ^ 2(A + 2 M ) £ ^_ 


Xfi 




2(A + 2/*) " 2(A + 2/ 


) *,0y 



(15) 



where we take the expression for e c a a from Eq. ([7]). We 
will look for linear combinations of the derivatives of the 
radial solutions of the bi-harmonic equation (fl~4"l) which 
are linear in the eigen-strain and go to zero at large dis- 
tance. In addition the terms must transform as a vector 
field. Let u c ' T and u c '° be the solutions to the Lame- 
Navicr equations. We then have: 



G« + A) 



(P + A) 



dX a dX~, 
d 2 u c ^ 



/'• 



/'•- 



dXpdXp 
d 2 u c '° 



dX a dX 1 ^dXpdXfs 



(16) 



Using the radial solutions of the Lame-Navier equa- 
tion we can construct the following combinations which 
transform as a vector and also go to zero as r — > oo. 



Using the identities 




<9 2 lnr n 

dXpdXp ~ ' 


d 2 (r 2 lnr) 

3x, d xr A ^+ A 


we see from Eq (|17[) 




a 2 ^ T 

dXpdXp 


T 9 3 lnr 

" A dx a dx v dx x 


and similarly 




^ , , 


, ^.*,r dHnr 



dX a dX~, 



dX a dX v 8X x 



(19) 



(20) 



(21) 



*,T 



<91nr 



ut T = te'ji ^^+Be 



*.T 



d 3 h: 



-Ce 



, T 9 3 (r 2 lnr) Plugging Eq J20), J2T) into Eq JIB), we get 



and 



Q ^ dX/j /37 dX a dXpdX 1 ^ 8X a dXpdX 1 

(17) 



t£ u = A e 



*,o<91nr / * >0 3 3 h 



-C 



t|0 d 3 (r 2 lnr) 



C 



A(X + m) 
4(A + 2^) 



Q ^ dx p ^ dx a dx p dx- t ^ dx a dx p dx- t 

(18) We can now rewrite Eq (flTj) as 
I 



(22) 



u c,T =Ae * : T^: +Bel: 



9 3 lr 



A(A + /x) *, T <9 3 (r 2 lnr) 



a ^ 0Jfy ^ BXadXpdX^ 4(A + 2fi) ^ dXadXpdX-y 

We now compute the following identities 

ainr _ X 



dXa 



a 3 h 



— 2r [X a 8p^ + XpSa^ + XySap) + 8A" Q A" / gX 7 
dX a dX j3 dX 1 ~ ^ 

<9 3 (r 2 lnr) _ 2r 2 (AT Q <5 / 3 7 + X p 5 ai + X^5 a p) - 4X a X p X^ 



dX a dXpdX 1 

Using these identities, we can now write down Eq (|23[) as 



L .r = ^ e ..J§_ 



"a/3 r 2 



2B A(X + n) 



and similarly 



i^'° = A e 



,o x P 

a/3 r 2 



r 4 2(X + 2n)r 2 



2B A{\ + v) 



*,T 



€ Pj ( X a3f3y + Xp5 al + X 1 5 a p) + 



8B A{\ + fi) 



e£ (X a 5p 1 + Xf;S ai + X 1 8 a p) + 



r 4 2(A + 2/x)r 2 _ 
Now using Eq. J3), we can rewrite equations J2S) and (j2"6")l as 

cT ((n + (k)fJ-AX a 

2(A + 2 M )r 2 



r 6 (A + 2/z)r 4 
8B' A'(A + /i) 



r 6 (A + 2fi)r 4 



tp-f X a XpX 1 



£p 1 X a XpX 1 



(23) 



(24) 



(25) 



(26) 



(27) 



u c '° = 



A fj, 



4B 



(A + 2fi)r 2 



*,o 



X 6 



8B 



A'(\ + n) \ , i0 



r 6 (A + 2/^)r 4 



Cp^XaXpXj 



(28) 



The complete solution for the displacement field in the matrix is then given as 



U a = Ua + < 



A p. 

2(X + 2fi)r 2 ' V(A + 2^)r 



45 



*,0 TU" 



T ) e a 0^P 



8B A(\ + /t)\ „ 



6 ' (A + 2 M )^i e ^ X ^^ 



(29) 



Now at r — a (the inclusion boundary), the form of solution (|29|l must match with the Eshelby solution (|15p. This 
implies, 



i4 = 



q 2 (A + /x) 



A = a 2 , B = - 



a 4 (A + /x) 
8(A + 2/i) 



Plugging equation ([50]) into Eq. (p?9"|) . we get: 
{X + ft) fa 2 ' 



a 2(X + 2n)\r 



(Cn + Cfc)ATa + 



2 V | ° 2 \,.,0j 



2 fi_iW°:*2^ 



r 2 y c /37 r 2 



Noting that 



*,o 



p X P 



(Cn " Cfc) 



(2n a (h ■ X) - X a ) 



(30) 



(31) 



(32) 



and 



e^X a X^X 7 = (C " Cfc) (2(n • X) 2 - r 2 )X a 



(33) 



we find that Eq. (|3ip finally becomes 



u c (X) 



(Cn -Cfe)(A + M) / a 



4(A + 2/i) V^ 2 

„2 



, (Cn + Cfc) y _ 

'(Cn-Cfc) "' VA + /i ' r 2 



— + %\ 



2 1-- 2(n-f) 2 -l X 



(34) 



where r = — . The Cartesian components of Eq. (|34[) used for visualizing the displacement field in space are given 
below: 



Clv ,_ (Cn~Ck)(X + ^) fa 2 
U ' AX) - 4(A + 2 M ) \? 



2 (Cn + Cfc^. _ 

(Cn-Cfc) V^ + M ' r 



2u, a 2 \. 

H — 2 ] (a;cos2<p + ysm20) 



a 2 \ /(a; 2 — y 2 )cos2(j> + 2xysin2(j) 



<(X) 



(C« -Cfc)(A + /i) /a 



4(A + 2/z) 



, (Cn+Cfc) 
'(Cn-Cfc) 



y- 



2a a , 

1 — - {xsn\2q> — ycos2<p) 

X + ji r 2 



2 1 



a 2 \ / (x 2 — y 2 )cos2(j} + 2xysin2<p 



(35) 



where the unit vector h makes an angle <fi with the x-axis. Taking derivatives of the displacement field 



du c a _ (C„ -Cfc)(A + /i) fa 
dX R 



-4 



4(A + 2j») 
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.X n Xf 



/'• 



A + ^j 



2(n • f)no, — 



(Cn " Cfc) 

X n \ Xfi 
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+ -o I 2n Q n^ - 5 Q/3 



2^)(2(n-r) 2 



Aq^g 



2 ) ( (w • r)hp - (h-r) 



r ) r 



2 1- 



2(n-rr-l <J a £ 



(36) 



and 



du % _ (Cn-Cfc)(A + M)/a 



dX a 4(A + 2^) \r 

4 



(Cn + Cfe) /r ^a^/3 



A* , a ~Wn/.~ ~\» -X^AX, 



A + , + ^> 
2u a 2 

A + /1 r z 

\2"*a \ -A/3 



^)( 2 n a n -5 a p)-4(l-2^)(2(n-rf-l)?4 l 



l -~2 ){{n-r)n a -{n-r) *-f \-*- + 2 M _ _ ] [ 2 (n ■ r)' - 



The constrained strain in the matrix can be written as 

1 ( du c a du% 



afi 2\dXp • dX a 
Using equations (|36|) and (j37|) . Eq. (f38|) becomes 



e*pW- 4(A + 2 M ) ^ 



2 ^" + Cfc J ^-2^ftl- 
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A + /i r^ 



(n • r) 



h a Xp npX 



_^_ + ^( M „ % _^)_ 4 ( 1 _ 2 ^( 2(<i .^_ 1 )^ 



4 1 



(n • r) 



hpX a n a Xp 



2(n-r) 
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2 1 



2(n-f) 2 -l )S a 



It is easy to see that the trace of € c al3 (X) is given as 



e nn\^> 



(Cn ~Ck)fifa 



(A + 2 M ) V" 



2(n-f) 2 - 1 



(37) 



(38) 



X n Xt- 



(39) 



(40) 



We can now calculate the constrained stress in the elastic matrix due to the deformed Eshelby inclusion. It follows 
from Hooke's law: 



a c ap (X) = 2^e c ap (X) + Xe^(X)S aP 
Plugging Eq. (|39[) in Eq. (|4Tj). we get the final expression for the constrained stress 

2 ^M( Sa ,-2^pL 



c f Y ^_ (Cn-Cfc)M(A + ^) /^Q 2 
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IV. ENERGY OF A/" ESHELBY INCLUSIONS 
EMBEDDED IN THE MATRIX 



The energy of the J\f Eshelby inclusions embedded in a where u a ^ = -7^ 



I 

where the superscript i denotes the index of the inclu- 
sion and m denotes the matrix. Eq (|43[) can be re-written 
using the definition of the strain e a p = (l/2)(u a ,p+up ta ), 



linear-elastic matrix m is given by the following expres- 
sion 

(43) 






E ^U/° 



/3 I "a,/3 "I" «/3,a 



dV 



v-I2iLi v$ 



<p[<,p + u^ a dV (44) 



Using the symmetry of the stress tensor, we obtain 

N 

7 af3 u P,a L 



i N r 1 /■ 



*=i J y 



v-XiLi v$ 



(45) 
If there are no body forces, we also have the identity 

<J a pUj3^ a = (a a /3Up) t a ~ O a 0,aUp = (a a pU/3),a (46) 

Thus we can write Eq. (|4"5j) as 

N 



HgX,M,/" + U 



^-/iii v o 



'a/3 ",9 



where e* 'a is the eigen-strain of the i th Eshelby inclu- 
sion and so on. Note that in the expression for the strain 
in the inclusion given by Eq, (f5"Tj) , we have subtracted the 
contribution of eigen-strain from the constrained strain 
leaving only the elastic contribution to calculate correctly 
the elastic contribution to the energy. Using these ex- 
pressions, the elastic energy of the system can be written 
from Eq. (|4T)|) as follows: 



.V 



(47) 

Using Gauss's theorem, we convert these volume inte- 
grals into area integrals to obtain 

^ = 2 E / , <^< dS ~ E 2 / , <^< dS 

A T " On A 1 "On 



i=l Js o ^ 



pup - <T a pUp )n a 

(52) 

Since the traction force has to be continuous at the 
inclusion boundary (Newton's III rd law), we have 



8=1 ' o 






(48) 



where n 1 and n°° are unit vectors both pointing out- 
wards from the surfaces bounding the inclusion volume 
Vq and the matrix boundary respectively. Eq. (|48)) can 
be re-written as follows 



E=\f < p ufdS + 5 E /. (<?4 ~ °%}V?)<dS 

J s°° ^ i=1 J si \ J 

= 2 a <xP € Pl X "l U a dS+ o2^ t [ a <xf3 U - ° a 

= \<^f a v + \y, l(°U4 - <f>uj)w a ds 

8=1 "'•SfN ' 

Thus we can re- write Eq. (fi"2"j) as 



'a/3 "-a — °ap' l a 

which gives us from Eq. (|52[) , 

M 

8=1 J b i. 

We also have from equations ([50]) and (fSl 



u /3- u /3 =- e ^ X « 



(53) 



i i N r 

E = ^XpV + o E / . </J W " u ™)^ ( 54 ) 



(55) 



/3"/3 I "a 



i=i 
<,(X)=^<^(X)+a- 

8=1 



On plugging this expression into Eq. ()54[) wc finally 

it 

4=1 •'^O 

= ^S^V - \ £ $ /^ (<^) dV 

I 1 ^ /■ 

= ^%>t%V - - £ %\ / <pk a dV 

L L i=l -^o 

1 1 ^ ._ 



2" ap-ap ' n / , • u ~pa a/3 

j=l 



(56) 



(/ 



nrf 



w = E*S( jr ) 



ft 



a/3 



(50) 



i=l 



where e^'l(.X') denotes the constrained strain at X in 
the matrix due to the i th Eshelby inclusion. We also have 
for locations X inside the inclusions 



where g 1 „ = (1/Vq ) J vi <J l a pdV . Using the expression 
for a l a Q from equation Eq. (|51[) . wc obtain 






e^W=^e^(X)+e^(X)-e^ + e- 



OO 





j& 



< P (x) = j2 °*i w + ^ w - <^ + < 



&i 



C(X) =J2u c a 3 (X)+u c a \X) - e$X p +uZ(51) 



Eq. (|57p is a far field approximation which assumes 
that r'-' 3> a. As r w -> a, clearly the spatial integrals 

contributing to a°\ must be computed explicitly and 
cannot be replaced by the single distance r y between 
the centers of the Eshelby inclusions i and j. 



Using equations (|56|) and ([57|) , we can write down the 
final form of the elastic energy expression. 

1 1 M 1 N 

+ \ E #Wttf - 5 E 4'>o E «#(»■") 

i— 1 i— 1 j^i 

= £mat + -Boo + -Besh + Ei nc (58) 

where each component of energy is defined as 

rp _oo ^oo t/ 

fi mat — 2 CT a^ e /9Q! K 

1 N 

Eoo = - g^ E e /3« y o 

i=l 

1 ^ 
^esft ~ 2 U^olP ~ a a fs) e /3a V 

i=\ 

1 N 

^=-2E^T^^) ( 59 ) 

Here the eigen-strain e*'l and volume Vq associated 
with any i Eshelby inclusion are given as 

V*- — 

V °~ 2 

e: ,, = (^^ (2< ^_^ )+ (Cn+0)^ (6Q) 

Also for a 2D material being loaded under uni-axial 
strain with free boundaries along y, we can write the 
form of the global stress tensor as 







(61) 



By Hooke's law, we get the expression for applied 
global stress tensor 



Gap = 2 M e 2^ + ^ e ^a/3 

Taking trace of Eq. (|62[) , we find 

1 „ 1 



2(A + /i) ™ 2(A + /i) 
Plugging Eq. JB3J in Eq. dHJ) , we find 



(62) 



(63) 



(64) 



where 7 is the external strain. In the following, we 
discuss these components of elastic energy as shown in 
Eq. (J59J) in detail: 



Emat ■ It is the elastic energy that would be present in 
the strained matrix in the absence of inclusions. Plugging 



equation (|6lT) and (f64|) into Eq. (f59|) . we get the following 
expression 



E,, 



2ii{\ + ii) 1 2 V 



(65) 



.Boo : It is the contribution to the elastic energy caused 
by the Eshelby inclusions themselves. Note that this term 
can make a negative contribution with respect to E mat . 
Again plugging equations (f6"Tj) , (f64|) and (|60j) into Eq. 
(1551). we find 



N 



E n 



+ 



-27tqV(A + / i)7 ^ j (C„ - Cfc) 
(A + 2/«) 

(Cn+Ck) 



E 



-(2(4)* 



1) 
(66) 



To find the orientation of each Eshelby inclusion with 
respect to the principal stress direction, we need to min- 
imize Eoo with respect to 8, where 8 = cos _1 (n?|.) = 
sin" 1 (n'). We thus get 



d f(Cn -Cfc) 



dd 



2 cos 2 8- I) 



(Cn+Ck) 







or 



(67) 



Hence each Eshelby inclusion must be oriented along the 
principal stress direction (or 8 = 0). 

E es h- It is the self energy required to create the Es- 
helby inclusion and it is always positive. To calculate 
this term, we need the following quantities: 



C . I 

E 7* 



= C Q/ 3 7 50 7 5fei£ fc ; 

(A + M) 



2(A + 2a*) 
(A + /i) 



2(Cn + Cfc)(A + /i)<W 



(2<4 - 6 a0 ) 



(68) 



In deriving Eq. (f6"8")l . we have used equations (JTUJ), ([5]) 
and ([5]). The eigen stress for the i th Eshelby inclusion 
can be written [using Eq. Q] as 

0% = (A+A*)(Cn + Cfc)*a/3 + /*(Cn-C*)(2nJ,n^-* a /s) (69) 

Combining this and the expression for eigen-strain from 
Eq. ([60]), Eq. ([59]) becomes 



£, 



7ra 
"2" 



JV 



C,2 \ *,2 



E/ *,Z cA \ 



ira 2 Af[i{\ + A*) 
4(A + 2/Lt) 



{2(c„ + a) 2 + (Cn - a) 2 } (70) 



Einc- This term arises due to the interaction between 
Eshelby inclusions in the "far field approximation" . We 
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note from Eq. (|59|) 

1 N 



i=l 
,,2 



&i 



<»J> 



E 4'>^(^')+e;»(r«) (71) 



Recalling the value of eigen-strain (Eq. ([3])) and con- 
strained stress (Eq. (HH) due to an i th Eshelby inclusion, 
we write down Eq. (|7ip (after simplifying): 



<V> (A + 2/x) Vr«V 



(Ck-C) (n*T«r + (r^T«^-l + 



(a - a) 5 



-4 



A + /i r u ' 



—2 4(n* • nJ)(n* • r^)(n^ ■ r«) - 2(n 4 • r^') 2 - 2(n-?' • r'y') 2 + 1 



+ 2 



2/x a 2 



A + /i r u 
a 2 



2(n J • niy - 1 ) - 4( 1 - 2— ) [ 2(n l ■ r^f - 1 ) ( 2(n> ■ r^f - 1 



-r 161 1 — I I (n i ■ ni)(n l ■ r l i){ni ■ r l i) - (n l ■ r^) 2 (ni ■ r^f 

In deriving Eq. fj72|> . we have used the following identities: 



(72) 



t ap 



x$ x% 



r i3' 



(Cn-Cfc) l-2(n*-r«)' 



(n? -r l i) 



niX^ »JX«\ X^Xp 



f l 3 f l 3 I ~,ij 

(Cn - Cfe) ( 2(n i ■ nJ'Xn* • r^)(iV ■ r' l o) - (n i ■ r l J) 2 - (ni ■ r l i) 2 + 



r iif 
1 



e$(2ft£nj - <M = (C„ - 00 ( 2(n< • r^) 2 - 1 
*,i X %X% (Cn + Cfc) , (Cn-Cfc) 



°0 r *i 2 



2(n* -r^y - 1 



E a/3 



[rh> ■ r^ ) r-i— -f 



2 /„ ? «ii "\2 



2(nJ -r'J) 



2 *ggg 



- (n 4 ■r i 3y(nJ ■ r^) 

£* a g$al3 = (Cn + Cfe) 



2(Cn - Cfe) (»»* • n»)(n* • r«)(i»» ■ r«) 



(73) 



To calculate the shear band angle with respect to the 
principal direction of strain, we have to minimize Ei nc 
with respect to 9. Assuming all eigen directions to be the 
same, i.e taking n % = ni = n, (n ■ r* J ') 2 = cos 2 6* = \ and 
-Sjnj- — > (far field approximation) we find, by putting 
— E — 0- 



1 1 (Cn + Cfe) 

2 4 (C„ - Cfe) 



or 



6> = cos 



-1 /I l (Cn + Cfc) 
V 2 4 (Cn - Cfe) 



(74) 



It is very easy to see from Eq. ([71)) that the area pre- 
serving case, such as pure shear implies Cn = — Cfe leading 
to the angle exactly equal to 45°. As can be expected, 
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other loading conditions may result in different values of 
the angle. The two extreme cases occur for |Cn/Cfcl — ► 
and |Cn/Cfc| —> oo. The first case corresponds to an angle 
of 30° and the second - to an angle of 60°. Therefore 
all experimental observations should fall between these 
two extreme universal limits. Following Ref.|8(, we find 
that indeed all the experimental data presented there fall 
within our theoretical limits. We return now to our simu- 
lations shown in Fig. [TJto rationalize the angles observed. 
In order to understand the angles observed in our simu- 
lations we need to figure how different loading conditions 
affect the values of ( n and Cfc- We find that in the case 
of extension [see Fig ((4])] , the outward displacement sig- 
nificantly dominates the inward displacement, realizing a 
higher ratio of |Cn/Cfc| as compared to the case of com- 
pression [see Fig. (J5J]. 

We attribute this asymmetry to the steeper rise in the 
repulsive core as compared to the weakly attractive tail 
in any generic inter-particle interaction potentials. To es- 
timate the ratio \( n /(k\, we first calculate average length 
of both the incoming and the outgoing vectors in a small 
region around the core of the plastic event. The ratio of 
these lengths then determines |Cn/Cfe|- For compression, 
we find that |C«/Cfe| ~ 1-15 and for extension, we find 
ICn/Cfcl ~ 4.05. Plugging these values in Eq. ([71)1 . we 
find the shear band angle of 46° for the compression and 
54° - for the extension; both are in very good agreement 
with the angles observed in the simulations presented in 
Fig. (|T|). In the following section, we explain how the 
present atomistic theory can predict the yield strain un- 
der such loading conditions. 



PREDICTING YIELD STRAIN 



49.9 



-47.7 




(a) 






Z*? 




-52.0 



52.1 




In terms of our atomistic model, the global yield be- 
comes possible if the formation of infinitely many Es- 
helby inclusions is energetically favorable (in the ther- 
modynamic limit). In other words, the system should be 
able to create a density p = Af/L of such inclusions where 
L is the global linear scale of the system. Assuming all 
eigen strains and the eigen directions to be the same, we 
have from equations (|66[) and (|70| . 



FIG. 4. (a) A fundamental plastic event during an AQS uniax- 
ial extension simulation of a 2-dimensional amorphous solid. 
Shown is the non-affine displacement field in the whole sys- 
tem, (b) a small region around the core of the event shown 
in (a). Again for clarity we show only the incoming and out- 
going arrows. The ratio 4 s1 ~ 4.05 is determined by the ratio 
of the average lengths of the outgoing and incoming arrows. 



-27ra 2 /i(A + n)^ n 

E °° = (IT^) N ' 



(75) 



and 



E PX h — 



7ra 2 /i(A + (i) 
4(A + 2/i) 



{2(C„ + C fe ) 2 + (C„-Cfc) 2 }Af. (76) 



Since the inclusions are localized in a strip of dimensions 
La, the energy density of these two terms is computed as 



•Boo + E esh -2irap(\ + ^hCn 



La (A + 2//) 

ira 2 p{\ + p) 



(77) 



4(A + 2 M ) 



{2(c»+ar + (Cn-a)> 



It is easy to check that these two terms are the only ones 
that are linear in the density p (other terms are either 
of order p° or p 2 ). Now for 7 < 7 Y this energy density 
increases with p. The only solution that minimizes the 
energy is the single inclusion with p = 0. The condition 
that identifies 7 Y requires the derivative of this energy 
density with respect to p to vanish. In other words, the 
coefficient of p in Eq. ([77]) should vanish. Solving for the 
value of 7 that satisfies this condition we find 



Cn 

1Y = ~ 






2V Cn 



(78) 



In order to determine the yield strain we need in ad- 
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FIG. 5. (a) A fundamental plastic event during an AQS uni- 
axial compression simulation of a 2-dimensional amorphous 
solid. Shown is the non-affine displacement field in the whole 
system, (b) a small region around the core of the event shown 
in (a). For clarity we show only the incoming and outgoing 
arrows. The ratio 4 s - « 1.15 is determined by the ratio of the 
average lengths of the incoming and outgoing arrows. 



dition to the ratio Cn/Cfe a l so the value of ( n . To com- 
pute the latter value, we find the best fit for the analytic 
expression of the clastic field produced by the Eshclby 
quadrupolar structure (Eq. (|34[0 to our numerical find- 
ings. Such a fitting procedure yields the values of (a, £„) 
as (0.9, -0.14) and (0.9, 0.09) for compression and exten- 
sion respectively. In Fig ([5]), we show our fits for both 
cases. 

Using Eq. (|75|) and the estimates of £„ and Cn/Cfci we 
estimate the yield strain for the case of compression and 



FIG. 6. (a) The fundamental plastic event of Fig. [Da observed 
during extension is modeled by an Eshelby inclusion. The 
best fit parameters ia found to be f n ~ 0.09, a — 0.9. (b) The 
same for the plastic event of Fig. [5b with a best fit given by 
C„ ~ -0.14, a = 0.9. 



extension as 






l7vh 


s6% , 


for compression 


l7vh 


s3% , 


for extension 



(79) 



This should be compared with the observed values of 
5.5% and 3.5% respectively in Fig. ^j. We consider the 
agreement quite satisfactory. 



VI. CONCLUSIONS 

Experimental observations and numerical simulations 
show that plastic phenomena in amorphous solids 
demonstrate essential asymmetry between the cases of 
uniaxial compression and extension. These asymmetries 



13 



are manifested in different angles the shear bands form 
with respect to the principal direction of stress and in 
very different values of the yield strain. The results pre- 
sented in this paper demonstrate that both asymmetries 
can be quantitatively described on the basis of atomistic 
theory of plastic events. We also derive analytically that 
the values of the shear band angles lie between 30° — 60° 
in good agreement with available experimental data. We 
reiterate the essential steps: one calculated the energy 
associated with M Eshelby inclusions, and in view of our 
athermal conditions minimizes this energy to find the se- 
lected distribution of inclusions. We find that for 7 < 7y 
the only solution that minimizes the energy is that 
containing a single Eshelby inclusion. At 7 = 7 < 7Y 
a new branch of solutions can open up, allowing for a 



density of inclusions to establish itself. The minimum 
energy is realized by a line of equi-distant inclusions 
that aligns with and angle 9 to the principal stress axis. 
The angle depends on the loading conditions as encoded 
by the eigenvalues of the Eshelby quadrupole. Only for 
simple shear we expect this angle to be 45°, while in 
general it is limited between 30° and 60°. Finally we 
computed analytically the yield-strain asymmetry under 
uniaxial loading conditions. A natural extension of our 
present work is a 3d generalization. Other possible 
directions would be to include finite temperature and 
strain rate effects. 
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